Ondrej Cernotik
23 February 2015
The evolution of a qubit and a thermal cavity coupled via a Jaynes-Cummings interaction is given by the stochastic master equation $$ \mathrm{d}\rho = -i[g(a\sigma_+ e^{i(\omega t+\phi)}+a^\dagger\sigma_- e^{-i(\omega t+\phi)})+\Delta a^\dagger a,\rho]\mathrm{d}t +\kappa(\bar{n}+1)\mathcal{D}[a]\rho\mathrm{d}t+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho\mathrm{d}t +\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)ae^{-i\delta t}-\bar{n}a^\dagger e^{i\delta t}]\rho\mathrm{d}W. $$ The equation is written in the interaction picture with respect to the qubit Hamiltonian $H_q = \frac{\omega}{2}\sigma_z$ and we use a detuning in the local oscillator to measure the signal appearing in the sideband of the output light.
Starting from this equation and choosing $\Delta = \omega = -\delta$, the Gaussian adiabatic elimination method [1] gives the effective dynamics $$ \mathrm{d}\rho_q = \frac{4g^2}{\kappa}\{(\bar{n}+1)\mathcal{D}[\sigma_-]+\bar{n}\mathcal{D}[\sigma_+]\}\rho_q\mathrm{d}t+ \frac{2g}{\sqrt{\kappa(2\bar{n}+1)}}\mathcal{H}[(\bar{n}+1)\sigma_-e^{-i(\phi+\pi/2)}-\bar{n}\sigma_+e^{i(\phi+\pi/2)}]\rho_q\mathrm{d}W. $$ With density operator expansion [2] for the measurement and projection operator method [3] for the deterministic part, we get $$ \mathrm{d}\rho_q = \frac{4g^2}{\kappa}\{(\bar{n}+1)\mathcal{D}[\sigma_-]+\bar{n}\mathcal{D}[\sigma_+]\}\rho_q\mathrm{d}t+ \frac{2g}{\sqrt{\kappa}}\mathcal{H}[\sigma_-e^{-i(\phi+\pi/2)}]\rho_q\mathrm{d}W. $$ See Ref. [1] for details on the derivation of these equations.
Since the full equation has time-dependent measurement operator and interaction Hamiltonian, we move to the rotating frame of the cavity field to obtain the time-independent equation $$ \mathrm{d}\rho = -ig[a\sigma_+ e^{i\phi}+a^\dagger\sigma_- e^{-i\phi},\rho]\mathrm{d}t +\kappa(\bar{n}+1)\mathcal{D}[a]\rho\mathrm{d}t+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho\mathrm{d}t +\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)a-\bar{n}a^\dagger]\rho\mathrm{d}W. $$ This equation is taken as the full solution with which the effective equations will be compared. The difference between the solutions is again given by the trace distance of the corresponding density operators (see the file Dispersive qubit readout.ipynb for details).
References
[1] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.
[2] A. Doherty and K. Jacobs, Physical Review A 60, 2700 (1999).
[3] C. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004).
In [1]:
%matplotlib inline
In [2]:
from qutip import qeye, destroy, sigmax, sigmay, sigmaz, sigmam, sigmap, tensor, smesolve, ket2dm, basis, thermal_dm
from matplotlib.pyplot import subplots, plot
from numpy import linspace, exp, sqrt, pi
In [3]:
def trace_dist(g, phase, kappa, nbar, Ntraj, NFock, Nstep, Nsub, tmax, rhoq0):
#***** function trace_dist *****
#This function calculates the trace distance between the exact dynamics
#and the approximate solutions (Gaussian adiabatic elimination, density
#operator expansion) for the Jaynes-Cummings interaction.
#Input parameters:
# g: qubit-cavity coupling,
# phase: interaction phase,
# kappa: intrinsic measurement (and decay) rate,
# nbar: thermal occupation,
# Ntraj: number of generated trajectories,
# NFock: Fock number cutoff,
# Nstep: number of steps for evolution,
# Nsub: number of substeps,
# tmax: maximum time,
# rhoq0: initial qubit state.
#Returns:
# tlist: list of times,
# D_gae: trace distance for Gaussian adiabatic elimination,
# D_doe: trace distance for density operator expansion.
#***** System operators *****
tlist = linspace(0, tmax, Nstep)
Id2 = qeye(2)
IdN = qeye(NFock)
a = tensor(Id2, destroy(NFock))
x = (a+a.dag())/sqrt(2.)
p = -1j*(a-a.dag())/sqrt(2.)
sx = tensor(sigmax(), IdN)
sy = tensor(sigmay(), IdN)
sz = tensor(sigmaz(), IdN)
sm = tensor(sigmam(), IdN)
sp = tensor(sigmap(), IdN)
#***** Full dynamics *****
H_full = g*(a*sp*exp(1j*phase)+a.dag()*sm*exp(-1j*phase))
c_op_full = [sqrt(2*kappa*nbar*(nbar+1)/(2*nbar+1))*x]
sc_op_full = [sqrt(kappa/(2*nbar+1))*((nbar+1)*a-nbar*a.dag())]
e_op_full = [sx, sy, sz]
rho0 = tensor(rhoq0, thermal_dm(NFock, nbar))
#***** Gaussian adiabatic elimination *****
H_gae = 0*Id2
c_op_gae = [(2*g/kappa)*sqrt(kappa*nbar*(nbar+1)/(2*nbar+1))*(sigmam()*exp(-1j*(phi+pi/2))+sigmap()*exp(1j*(phi+pi/2)))]
sc_op_gae = [(2*g/sqrt(kappa*(2*nbar+1)))*((nbar+1)*exp(-1j*(phi+pi/2))*sigmam()-nbar*exp(1j*(phi+pi/2))*sigmap())]
e_op_gae = [sigmax(), sigmay(), sigmaz()]
#***** Density operator expansion *****
H_doe = H_gae
c_op_doe = [(2*g/sqrt(kappa))*nbar*sigmam(), (2*g/sqrt(kappa))*nbar*sigmap()]
sc_op_doe = [(2*g/sqrt(kappa))*exp(-1j*(phi+pi/2))*sigmam()]
e_op_doe = e_op_gae
#***** Generating trajectories *****
D_gae = 0
D_doe = 0
for i in xrange(0,Ntraj):
print i
sol_full = smesolve(H_full, rho0, tlist, c_op_full, sc_op_full, e_op_full, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', store_measurement=True)
sol_gae = smesolve(H_gae, rhoq0, tlist, c_op_gae, sc_op_gae, e_op_gae, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', noise=sol_full.noise)
sol_doe = smesolve(H_doe, rhoq0, tlist, c_op_doe, sc_op_doe, e_op_doe, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', noise=sol_full.noise)
s1x = sol_full.expect[0]
s1y = sol_full.expect[1]
s1z = sol_full.expect[2]
s2x = sol_gae.expect[0]
s2y = sol_gae.expect[1]
s2z = sol_gae.expect[2]
s3x = sol_doe.expect[0]
s3y = sol_doe.expect[1]
s3z = sol_doe.expect[2]
D_gae = D_gae+sqrt((s1x-s2x)**2+(s1y-s2y)**2+(s1z-s2z)**2)/2.
D_doe = D_doe+sqrt((s1x-s3x)**2+(s1y-s3y)**2+(s1z-s3z)**2)/2.
D_gae = D_gae/Ntraj
D_doe = D_doe/Ntraj
return tlist, D_gae, D_doe
In [4]:
g = 0.1
phi = 0.
kappa = 1.
nbar = 0.5
Ntraj = 20
NFock = 8
Nsteps = 1000
Nsub = 250
tmax = 50
t, D_gae, D_doe = trace_dist(g, phi, kappa, nbar, Ntraj, NFock, Nsteps, Nsub, tmax, ket2dm((basis(2,0)+basis(2,1)).unit()))
In [5]:
fig, axes = subplots(figsize=(12,8))
axes.plot(t, D_gae, label='Gaussian ad. elimination')
axes.plot(t, D_doe, label='Density op. expansion')
axes.set_xlabel('Time')
axes.set_ylabel('Trace distance')
axes.legend(loc=0)
Out[5]:
In [6]:
D_gae_time = D_gae.sum()/Nsteps
print D_gae_time
In [7]:
D_doe_time = D_doe.sum()/Nsteps
print D_doe_time